clc; clear;
cd 'data'

lambda_data = csvread('lambda.csv');
rho_data = csvread('rho.csv');
Y_data = csvread('Y.csv');

cd ..

I = 44;

alpha = 1; 
theta = 2.5; sigma = 3;
YEAR = 15;

for t = 1:YEAR
   lambda{t} = reshape(lambda_data(:,t),I,I);    
   rho{t} = 1-reshape(rho_data(:,t),I,I);    
   Y(:,t) = Y_data(:,t);
   rho_autarky(:,t) = max(0,1-(1-diag(rho{t}))./diag(lambda{t}));
end

for t = 1:YEAR
   LL(:,t) = diag(lambda{t});    
end

Real_change  = LL.^(-1/theta);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for t = 1:YEAR
 
Gini(:,t) = alpha.*(sigma-1)./(sigma*theta) - 2.*alpha.*(theta-(sigma-1))./(sigma*theta).*(sigma-1)./(4*theta-2*(sigma-1))...
    .* sum(lambda{t}.*(1-rho{t}).*repmat(Y(:,t)',I,1)./repmat(Y(:,t),1,I),2);
Gini_autarky(:,t) = alpha.*(sigma-1)./(sigma*theta) -  2.*alpha.*(theta-(sigma-1))./(sigma*theta).*(sigma-1)./(4*theta-2*(sigma-1)).*(1-rho_autarky(:,t));

s = 0.95; star = rho{t};

Lorenz5(:,t) = (1-alpha.*(sigma-1)./(sigma*theta)).*s + alpha.*(theta-(sigma-1))./(sigma*theta).*sum((ones(I,I).*(s>star))...
    .*lambda{t}.*repmat(Y(:,t)',I,1)./repmat(Y(:,t),1,I).*(theta./(theta-(sigma-1)).*(1-((1-s)./(1-rho{t})).^((1-sigma)/theta+1))...
    -1 + (1-s)./(1-rho{t})),2);

Share5(:,t) = 1 - Lorenz5(:,t);
 
s = 0.99; 
Lorenz1(:,t) = (1-alpha.*(sigma-1)./(sigma*theta)).*s + alpha.*(theta-(sigma-1))./(sigma*theta).*sum((ones(I,I).*(s>star))...
    .*lambda{t}.*repmat(Y(:,t)',I,1)./repmat(Y(:,t),1,I).*(theta./(theta-(sigma-1)).*(1-((1-s)./(1-rho{t})).^((1-sigma)/theta+1))...
    -1 + (1-s)./(1-rho{t})),2);
Share1(:,t) = 1 - Lorenz1(:,t);


s = 0.95; star_autarky = rho_autarky(:,t);
Lorenz5_autarky(:,t) = (1-alpha.*(sigma-1)./(sigma*theta)).*s + alpha.*(theta-(sigma-1))./(sigma*theta).*(s>star_autarky)...
   .*(theta./(theta-(sigma-1)).*(1-((1-s)./(1-rho_autarky(:,t))).^((1-sigma)/theta+1))...
    -1 + (1-s)./(1-rho_autarky(:,t)));
Share5_autarky(:,t) = 1 - Lorenz5_autarky(:,t);

s = 0.99; 
Lorenz1_autarky(:,t) = (1-alpha.*(sigma-1)./(sigma*theta)).*s + alpha.*(theta-(sigma-1))./(sigma*theta).*(s>star_autarky)...
   .*(theta./(theta-(sigma-1)).*(1-((1-s)./(1-rho_autarky(:,t))).^((1-sigma)/theta+1))...
    -1 + (1-s)./(1-rho_autarky(:,t)));
Share1_autarky(:,t) = 1 - Lorenz1_autarky(:,t);
end

Gini_change = 100*(-1+Gini./Gini_autarky);
S5_change = 100*(-1+Share5./Share5_autarky);
S1_change = 100*(-1+Share1./Share1_autarky);
R_change = 100*(-1 + Real_change);
 
Gini_change = Gini_change(:,end)-Gini_change(:,1);
S5_change =   S5_change(:,end)-S5_change(:,1);
S1_change =   S1_change(:,end)-S1_change(:,1);
R_change =    R_change(:,end)-R_change(:,1);

iso = {"AUS","AUT","BEL","BGR","BRA","CAN","CHE","CHN","CYP","CZE","DEU","DNK","ESP","EST","FIN","FRA","GBR","GRC","HRV","HUN","IDN","IND","IRL","ITA","JPN","KOR","LTU","LUX","LVA","MEX","MLT","NLD","NOR","POL","PRT","ROU","ROW","RUS","SVK","SVN","SWE","TUR","TWN","USA"};
plot(R_change, Gini_change,'w.'); hold on;
text(R_change, Gini_change,iso,'FontSize', 14,'Color', [0.6350 0.0780 0.1840]);
xlim([-2 14]);
ylim([-2 8]);
set(gca,'FontSize',20);
xlabel('Relative Change in Welfare (2000-2014)','FontSize',30, 'Interpreter','latex');
ylabel('Relative Change in Gini (2000-2014)','FontSize',30, 'Interpreter','latex');
grid on;

figure;

plot(R_change, S5_change,'w.'); hold on;
text(R_change, S5_change,iso,'FontSize', 14,'Color', [0 0.4470 0.7410]);
xlim([-2 14]);
ylim([-2 8]);
set(gca,'FontSize',20);
xlabel('Relative Change in Welfare (2000-2014)','FontSize',30, 'Interpreter','latex');
ylabel('Relative Change in Share of Top 5\%','FontSize',30, 'Interpreter','latex');
grid on;

%=========================================================================%

%% Numbers in the paper
mean(Gini_change, 'all')
mean(S5_change, 'all')
mean(S1_change, 'all')
mean(R_change, 'all')

